% smoothed conditional volatilities for cumulative inflation
clear

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
year = 1850:1:2013;
%select_years = [1864, 1904, 1921, 1959, 1980, 1994, 2009];
select_years = [15, 55, 72, 110, 131, 145, 160];
NY = size(select_years,2);
HZN = 10; % longest forecast horizon

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% catalog of data files
DFILE(1,:) = ['swuc_wpi_ar1_01'];
DFILE(2,:) = ['swuc_wpi_ar1_02'];
DFILE(3,:) = ['swuc_wpi_ar1_03'];
DFILE(4,:) = ['swuc_wpi_ar1_04'];
DFILE(5,:) = ['swuc_wpi_ar1_05'];
DFILE(6,:) = ['swuc_wpi_ar1_06'];
DFILE(7,:) = ['swuc_wpi_ar1_07'];
DFILE(8,:) = ['swuc_wpi_ar1_08'];
DFILE(9,:) = ['swuc_wpi_ar1_09'];
DFILE(10,:) = ['swuc_wpi_ar1_10'];
DFILE(11,:) = ['swuc_wpi_ar1_11'];
DFILE(12,:) = ['swuc_wpi_ar1_12'];
DFILE(13,:) = ['swuc_wpi_ar1_13'];
DFILE(14,:) = ['swuc_wpi_ar1_14'];
DFILE(15,:) = ['swuc_wpi_ar1_15'];
DFILE(16,:) = ['swuc_wpi_ar1_16'];
DFILE(17,:) = ['swuc_wpi_ar1_17'];
DFILE(18,:) = ['swuc_wpi_ar1_18'];
DFILE(19,:) = ['swuc_wpi_ar1_19'];
DFILE(20,:) = ['swuc_wpi_ar1_20'];
DFILE(21,:) = ['swuc_wpi_ar1_21'];
DFILE(22,:) = ['swuc_wpi_ar1_22'];
DFILE(23,:) = ['swuc_wpi_ar1_23'];
DFILE(24,:) = ['swuc_wpi_ar1_24'];
DFILE(25,:) = ['swuc_wpi_ar1_25'];
DFILE(26,:) = ['swuc_wpi_ar1_26'];
DFILE(27,:) = ['swuc_wpi_ar1_27'];
DFILE(28,:) = ['swuc_wpi_ar1_28'];
DFILE(29,:) = ['swuc_wpi_ar1_29'];
DFILE(30,:) = ['swuc_wpi_ar1_30'];
DFILE(31,:) = ['swuc_wpi_ar1_31'];
DFILE(32,:) = ['swuc_wpi_ar1_32'];
DFILE(33,:) = ['swuc_wpi_ar1_33'];
DFILE(34,:) = ['swuc_wpi_ar1_34'];
DFILE(35,:) = ['swuc_wpi_ar1_35'];
DFILE(36,:) = ['swuc_wpi_ar1_36'];
DFILE(37,:) = ['swuc_wpi_ar1_37'];
DFILE(38,:) = ['swuc_wpi_ar1_38'];
DFILE(39,:) = ['swuc_wpi_ar1_39'];
DFILE(40,:) = ['swuc_wpi_ar1_40'];

NF = size(DFILE,1);
NB = 21;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load posterior sample
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% variance of innovations to log volatilities
% load WA
load(DFILE(NB,:),'WD')
[P,N,N] = size(WD);
WA = zeros((NF-NB+1)*P,N,N);
WA(1:P,:,:) = WD;
clear WD 
for i = NB+1:NF,
  j = i - NB + 1; 
  load(DFILE(i,:),'WD')
  WA((j-1)*P+1:j*P,:,:) = WD;
  clear WD;
end
NMC = size(WA,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  compute RA and QA
load(DFILE(NB,:),'vD')
[N,T,P] = size(vD);
RA = zeros(T,(NF-NB+1)*P);
QA = zeros(T,(NF-NB+1)*P);
RA(:,1:P) = vD(1,:,:);
QA(:,1:P) = vD(2,:,:);
clear vD 
for i = NB+1:NF,
  j = i - NB + 1; 
  load(DFILE(i,:),'vD')
  RA(:,(j-1)*P+1:j*P) = vD(1,:,:);
  QA(:,(j-1)*P+1:j*P) = vD(2,:,:);
  clear vD;
end

% hiddent states
% compute SA
load(DFILE(NB,:),'SD')
[P,K,T] = size(SD);
SA1 = zeros((NF-NB+1)*P,T); % \pi
SA2 = zeros((NF-NB+1)*P,T); % \mu
SA1(1:P,:,:) = squeeze(SD(:,1,:));
SA2(1:P,:,:) = squeeze(SD(:,2,:));
clear SD 
for i = NB+1:NF,
  j = i - NB + 1; 
  load(DFILE(i,:),'SD')
  SA1((j-1)*P+1:j*P,:,:) = squeeze(SD(:,1,:));
  SA2((j-1)*P+1:j*P,:,:) = squeeze(SD(:,2,:));
  clear SD;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NMC = size(RA,2);
hgrid = 0:.01:1.5; % grid for histograms
HCSV = zeros(size(hgrid,2),2,T); % histograms for conditional smoothed volatilities
IQR1 = zeros(3,2,T); % conditional std
IQR2 = zeros(3,2,T); % conditional root mean square
CSV = zeros(2,NMC,NY);
CRMS = zeros(2,NMC,NY);

for tt = 1:T,
    
    tt
    c1m = zeros(HZN,NMC);
    c2m = zeros(HZN,NMC);
    sa2 = SA2(:,tt);
    ra = RA(tt+1,:);
    qa = QA(tt+1,:);
    matlabpool local 8
    parfor ii = 1:NMC,
          [c1m(:,ii),c2m(:,ii)] = conditional_2nd_moment_cum_inflation_sw(sa2(ii),ra(ii),qa(ii),squeeze(WA(ii,:,:)),HZN);
    end
    matlabpool close
    cstd = (c2m([5 10],:) - c1m([5 10],:).^2).^.5;
    Ncstd = hist(cstd',hgrid);
    HCSV(:,:,tt) = Ncstd/NMC;
    ss = sort(cstd');
    IQR1(:,:,tt) = ss([.25 .5 .75]*NMC,:);
    crms = c2m([5 10],:).^.5;
    ss = sort(crms');
    IQR2(:,:,tt) = ss([.25 .5 .75]*NMC,:);
    
    % store smoothed volatilities for selected years
    [mm,ind] = min(abs(tt*ones(1,NY) - select_years));
    if mm == 0,
        CSV(:,:,ind) = cstd;
        CRMS(:,:,ind) = c2m([5 10],:).^.5; 
    end
   
end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% posterior IQR for smoothed volatilities
figure
subplot(2,2,1)
plot(year,squeeze(IQR1(2,1,:)),'-k','Linewidth',2); hold on;
plot(year,squeeze(IQR1([1 3],1,:)),':k','Linewidth',2); hold off;
legend('Median','Interquartile range','Location','Northeast')
title('5 years ahead','Fontsize',16)
set(gca,'Fontsize',16)
axis([1850 2013 0 0.8])
grid

subplot(2,2,2)
plot(year,squeeze(IQR1(2,2,:)),'-k','Linewidth',2); hold on;
plot(year,squeeze(IQR1([1 3],2,:)),':k','Linewidth',2); hold off;
title('10 years ahead','Fontsize',16)
set(gca,'Fontsize',16)
axis([1850 2013 0 0.8])
grid

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% posterior IQR for crms 
figure
subplot(2,2,1)
plot(year,squeeze(IQR2(2,1,:)),'-k','Linewidth',2); hold on;
plot(year,squeeze(IQR2([1 3],1,:)),':k','Linewidth',2); hold off;
legend('Median','Interquartile range','Location','Northeast')
title('5 years ahead','Fontsize',16)
set(gca,'Fontsize',16)
axis([1850 2013 0 1.1])
grid

subplot(2,2,2)
plot(year,squeeze(IQR2(2,2,:)),'-k','Linewidth',2); hold on;
plot(year,squeeze(IQR2([1 3],2,:)),':k','Linewidth',2); hold off;
title('10 years ahead','Fontsize',16)
set(gca,'Fontsize',16)
axis([1850 2013 0 1.1])
grid


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CSV before and after WWII 
ECSV = squeeze(mean(CSV,2)); % mean smoothed volatilities
rmsv5 = ECSV(1,:)/ECSV(1,1); % 1864 is base year
rmsv10 = ECSV(2,:)/ECSV(2,1);
[rmsv5(4:7)' rmsv10(4:7)']

rmsv5 = ECSV(1,:)/ECSV(1,2); % 1903 is base year
rmsv10 = ECSV(2,:)/ECSV(2,2);
[rmsv5(4:7)' rmsv10(4:7)']

rmsv5 = ECSV(1,:)/ECSV(1,3); % 1921 is base year
rmsv10 = ECSV(2,:)/ECSV(2,3);
[rmsv5(4:7)' rmsv10(4:7)']

% PROBABILITY OF A DECLINE RELATIVE TO A BASE YEAR
CSV5 = squeeze(CSV(1,:,:));
CSV10 = squeeze(CSV(2,:,:));
PR = zeros(4,2);
for yy = 4:7,
    rcstd = CSV5(:,yy)./CSV5(:,1); % 1864
    DD5 = zeros(size(rcstd));
    DD5(rcstd<1) = 1;
    PR(yy-3,1) = mean(DD5');

    rcstd = CSV10(:,yy)./CSV10(:,1);
    DD10 = zeros(size(rcstd));
    DD10(rcstd<1) = 1;
    PR(yy-3,2) = mean(DD10');
end
PR

PR = zeros(4,2);
for yy = 4:7,
    rcstd = CSV5(:,yy)./CSV5(:,2); % 1903
    DD5 = zeros(size(rcstd));
    DD5(rcstd<1) = 1;
    PR(yy-3,1) = mean(DD5');

    rcstd = CSV10(:,yy)./CSV10(:,2);
    DD10 = zeros(size(rcstd));
    DD10(rcstd<1) = 1;
    PR(yy-3,2) = mean(DD10');
end
PR
 
PR = zeros(4,2);
for yy = 4:7,
    rcstd = CSV5(:,yy)./CSV5(:,3); % 1921
    DD5 = zeros(size(rcstd));
    DD5(rcstd<1) = 1;
    PR(yy-3,1) = mean(DD5');

    rcstd = CSV10(:,yy)./CSV10(:,3);
    DD10 = zeros(size(rcstd));
    DD10(rcstd<1) = 1;
    PR(yy-3,2) = mean(DD10');
end
PR

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% comparing post WWII with pre WWI: CRMS
ECRMS = squeeze(mean(CRMS,2)); % mean smoothed RMS
rm5 = ECRMS(1,:)/ECRMS(1,1); % 1864 is base year
rm10 = ECRMS(2,:)/ECRMS(2,1);
[rm5(4:7)' rm10(4:7)']

rm5 = ECRMS(1,:)/ECRMS(1,2); % 1903 is base year
rm10 = ECRMS(2,:)/ECRMS(2,2);
[rm5(4:7)' rm10(4:7)']

rm5 = ECRMS(1,:)/ECRMS(1,3); % 1921 is base year
rm10 = ECRMS(2,:)/ECRMS(2,3);
[rm5(4:7)' rm10(4:7)']

% PROBABILITY OF A DECLINE RELATIVE TO A BASE YEAR
CRMS5 = squeeze(CRMS(1,:,:));
CRMS10 = squeeze(CRMS(2,:,:));
PR = zeros(4,2);
for yy = 4:7,
    rcrms = CRMS5(:,yy)./CRMS5(:,1); 
    DD5 = zeros(size(rcrms));
    DD5(rcrms<1) = 1;
    PR(yy-3,1) = mean(DD5');

    rcrms = CRMS10(:,yy)./CRMS10(:,1);
    DD10 = zeros(size(rcrms));
    DD10(rcrms<1) = 1;
    PR(yy-3,2) = mean(DD10');
end
PR

PR = zeros(4,2);
for yy = 4:7,
    rcrms = CRMS5(:,yy)./CRMS5(:,2); 
    DD5 = zeros(size(rcrms));
    DD5(rcrms<1) = 1;
    PR(yy-3,1) = mean(DD5');

    rcrms = CRMS10(:,yy)./CRMS10(:,2);
    DD10 = zeros(size(rcrms));
    DD10(rcrms<1) = 1;
    PR(yy-3,2) = mean(DD10');
end
PR

PR = zeros(4,2);
for yy = 4:7,
    rcrms = CRMS5(:,yy)./CRMS5(:,3); 
    DD5 = zeros(size(rcrms));
    DD5(rcrms<1) = 1;
    PR(yy-3,1) = mean(DD5');

    rcrms = CRMS10(:,yy)./CRMS10(:,3);
    DD10 = zeros(size(rcrms));
    DD10(rcrms<1) = 1;
    PR(yy-3,2) = mean(DD10');
end
PR
